{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc4 as pm\n",
    "import tensorflow as tf\n",
    "from tensorflow_probability import bijectors\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Logp calculation for linear regression \n",
    "@pm.model(auto_name=True)\n",
    "def linreg(n_points=100):    \n",
    "    # Define priors\n",
    "    sigma = pm.HalfNormal(sigma=1., bijector=bijectors.Identity())\n",
    "    intercept = pm.Normal(mu=0, sigma=.1)\n",
    "    x_coeff = pm.Normal(mu=0, sigma=1.)\n",
    "    x = np.linspace(-5, 5, n_points)\n",
    "    \n",
    "    # Define likelihood\n",
    "    y = pm.Normal(mu=intercept + x_coeff * x, sigma=sigma)\n",
    "    \n",
    "model = linreg.configure()\n",
    "\n",
    "forward_sample = model.forward_sample()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0xb2d050f60>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGOFJREFUeJzt3W+MHVd5x/HfE7M0S/mzVDFCWSfYSMVAaxTTJaJatRQnkFCnieUihSIQLS+sIkCJBaYb8qZ9U29riVAJ1MpCVEhEjRGkbkugJpHTSrVEyhobQkhcpQSINyAWwQJqVomdPH2xe5Pr65l7Z+6cmXNm7vcjIbHr65kzRjxz7nOe8xxzdwEAuuOS2AMAAIRFYAeAjiGwA0DHENgBoGMI7ADQMQR2AOgYAjsAdAyBHQA6hsAOAB3zghg3veyyy3zr1q0xbg0ArXXy5MmfuvvmUZ+LEti3bt2qpaWlGLcGgNYysx8U+RypGADoGAI7AHQMgR0AOobADgAdEySwm9mMmX3RzB4xs4fN7HdDXBcAUF6oqpi/k/Tv7v5OM3uhpBcFum4tjp5a1qFjZ/TE6poun5nWgeu2a8/O2djDAoAgKgd2M3uppN+X9KeS5O5PS3q66nXLKBOoj55a1m13P6i1c89IkpZX13Tb3Q9KUrDgzosDQEwhUjGvlrQi6R/N7JSZfcbMfj3AdQvpBerl1TW5ng/UR08tZ37+0LEzzwX1nrVzz+jQsTNRxgMAoYUI7C+Q9EZJf+/uOyX9n6SFwQ+Z2T4zWzKzpZWVlQC3XVc2UD+xulbq93WPBwBCCxHYz0o66+4PbPz8Ra0H+gu4+2F3n3P3uc2bR+6ILaxooD56alnzi8eVd3T35TPTjY4HAOpSObC7+48lPW5m2zd+dY2k71a9blF5Abn/9/3pkSzTU5t04LrtmX9Wx3gAoE6h6tg/LOlOM/u2pKsk/XWg64504Lrtmp7adMHvBgN1VnqkZ3ZmWgf37gi2uFlkPABQpyDlju5+WtJciGuV1QvIw6pQ8tIgJunEwq7GxwMAdYrS3TG0PTtnhwbOy2emM9MwdaVHRo0HAOo0ES0FSI8AmCSdmLGPMk56pIlNRmxkAlCHiQjsUrn0SFO7U+u+B4DJNBGpmLKa2GTERiYAdSGwZ2hikxEbmQDUhcCeoYlNRmxkAlAXAnuGJqpoqNQBUJeJWTwto4lNRmxkAlAXc89ri1Wfubk5X1paavy+ANBmZnbS3Ufu8mfGnjhq3QGURWBPGLXuAMbB4mnCqHUHMA4Ce8KodQcwDgJ7wqh1BzAOAnvCqHUHMA4WTxNGrTuAcRDYE8ehHQDKIhUDAB1DYAeAjiGwA0DHENgBoGMI7ADQMcGqYsxsk6QlScvufkOo66aMBl0AUhSy3PEWSQ9LemnAayaLBl0AUhUkFWNmWyTtlvSZENdrAxp0AUhVqBn7JyV9TNJLAl2vMUXSKVmfoUEXgFRVnrGb2Q2SfuLuJ0d8bp+ZLZnZ0srKStXbBtFLpyyvrsn1fDrl6KnlkZ+ZedFU5jVp0AUgthCpmHlJN5rZ9yXdJWmXmX1+8EPuftjd59x9bvPmzQFuW12RdEreZ9xFgy4ASaoc2N39Nnff4u5bJb1L0nF3f0/lkTWgSDol7zO/WDung3t3aHZmWiZpdmZaB/fuYOEUQHStbQIWotTw8plpLWcEbpc0v3hcB67bnvuZy2emJ6pBF6WdQHsE3aDk7v/RRA17kdx4EVn9znt613zrazdPfMol1L83gGa0cudpqFLDPTtnn0unZFk794zuf2Rl4lMulHYC7dLKVEzIUsNeOmXbwj3ynGv2p1x6KYn9R05PTEqC0k6gXVo5Y6/jLNAi15zUlARnrwLt0srAPuws0KOnljW/eFzbFu7R/OLxwkG3yPmiqaYkxn3mojh7FWiXVqZi8s4ClTR2/5Yi54ummJJoomcNZ68C7WLuWZnles3NzfnS0lLw684vHs8sTZydmdaJhV3JX38cIcdESSOQNjM76e5zoz7XmlRMkXRD3TPqFFMSoZ55UtcPgC5qRWAvGnTqXuTrL48sWvpYd/471DOnun4AoLxW5NiHBZ3+oHrguu0X5Jul8DPqvN2mWWkMafycf5ase4R65hTXDwCMpxUz9qJBZ3BGPTM9pUunLtH+I6drmS335H2j+Kt/eyjYLDjvHpKCbKCipBHojlbM2If1axnUm1E3ecJR3jeKwd/1jDMLHvat5cTCrsrP1MS3HQDNaMWMfZxFyyZzxmUD9Tiz4LpTJeOsHwBIUytm7OPUUTeZM877RjEzPaWnzj8bZBZc5lvLuKp2q6RcEkhDKwK7VD7oNBEIe/LSGH95429JCrOxJ+sepvUUU6/FcMwgyuHeQDpaE9jLajJnPOobRYjA1n+P5dU1mfRc07Kmg2jWzLxo5VIT+OaASdepnaeDuvp/8Jg7YAdn5tL6CzNvodgkPba4u9Yx9csbH+sF6IKiO087O2OXqueMUxWz5jxvZr7JTM9kTBKaLpdM6ZsDEEsrqmJwoZg153kvj2fck2i3wEYrgMCejDKtB+poW1xU3sujVx4Zu1ySjVZAx1MxbVG2oqSOtsVFDVuUTiH1xUYroOOLp20RajG0qUXV1BelUx8fMC4WT1skVF646HWqBr4UZubDpD4+oG4E9oZlBdVQm6mKXKdM2oeZL9BOlRdPzewKM7vfzB42s4fM7JYQA+uivA6Nb33t5iAVJSHPbeXgDaC9QlTFnJf0EXd/naQ3S/qgmb0+wHU7Jy+o3v/ISpCKkiKNvIqmazh4A2ivyqkYd/+RpB9t/PdfmdnDkmYlfbfqtbtmWFANlRcedZ2iaR/qwYH2ClrHbmZbJe2U9EDI63ZFCjXWRVsgpzBWAOMJFtjN7MWSviTpVnf/Zcaf7zOzJTNbWllZCXXbVknhMOyip0ylMFYA4wlSx25mU5K+LOmYu39i1OcnuY49pUqTUQ2zUhorgOJ17JUDu5mZpM9J+pm731rk70xyYE9JzC6RAMprcoPSvKT3SnrQzE5v/O7j7v6VANdGjWIvkPKNAKhHiKqY/9J62220TIiNUeMGZ05cAupDd8cJVnWBtMompqbq5OvudgmkiMA+wYpsaBqmSnBuIg3E7llMKnrFTLgqG6OqBOcmDhvnNCVMKmbsGFuVTUxN1MnHXhwuglQR6kBgx9iqBOeqaaAiUt89S6oIdSEVk6C2lAHmneRUdKx1901P/TQlUkWoC4E9MW0rAywSnGO9qKq+eOrWhlQR2onAnpiuzeLqelEVfVnkvXhS+FbUxAIyJhM59sR0bRZXR716Vm56/5HT2lpwAbKJ3HaRRVEaraEuBPbEpL7gV1YdL6qsl0Wv41GRIF335qiiL44mFpAxmUjFJCb1Bb+ieqmOvBZzw15Uo9Iko14Ko1JXdX8rKpNO4+Bt1IHAnpjUF/yKyGoH3G/Yi6pITj4vN91vWJCuO7fdtXQa2odUTIL27JzViYVdemxxt04s7GpVUJeyZ6w9eQd7DPu7g2mSrNz0oGFBuu7cdtfSaWgfAnuHxdrVOGxm+tT5Z/XzJ8/l5p6LzHb7c9PSxa1FRwXpunPbLIoiNlIxHRWzHj4v1bHJbGTuuWiapD83PU7pYp257S6k09BuQY7GK4sTlOoX83SkvCP38tIzJumxxd1D/y7VIkCzJyghQTEX8PJmrIeOnRk5G2e22w4pbPBCPgJ7R8Xe1ZiX6ihSykkJYNra1vZiErF42lFZC3im9f8TxmoPy4acbmjq9CuMjxl7R/WnNJZX12S6eHdm/+eaHBeBvN2o008fM/YO69XDz85MX7QDlBkWxkWdfvoI7BOAGRZCok4/faRiJkDshdQ2odpjNCqX0kdgnwBdaSxWN6o9imOtJG1BUjFmdr2ZnTGzR81sIcQ1EU4bqlFSONSZag90ReUZu5ltkvRpSW+TdFbSN8zsX939u1WvjXBSnmGlMlNmLQJdEWLGfrWkR939e+7+tKS7JN0U4LqYEKnMlKn2QFeECOyzkh7v+/nsxu8uYGb7zGzJzJZWVlYC3BZdkcpMmWoPdEWIwD7YNVXSxQfnuPthd59z97nNmzcHuC26IpWZchvWIoAiQlTFnJV0Rd/PWyQ9EeC6mBApVe2kvBYBFBUisH9D0m+a2TZJy5LeJendAa6LCUFddHXU36Nf5cDu7ufN7EOSjknaJOmz7v5Q5ZFhoqQ+U045cKZSVYR0BNmg5O5fkfSVENcCUtEL5ik1UcsyrKoohfGhefSKATL0ZsG9VgwpN1FLpaoI6aClAJAhaxY8KJXAOawXUMopJNSHGTuQoUjQTmXj0rBDVfYfOa3l1TW5nk8hxWjXgGYR2IEMo4J2ShuX+uvvJV2wHpByCgn1IbBjYpRpNJY3C5bS3Lg07FCVQamkkFAfcuyYCGVLAkPW1teV5866bptSSKiPuY96v4c3NzfnS0tLjd8X7VYlQM4vHs9cYJydmdaJhV2hh/qcwReKtJ7GyZrxl3m+vOteOnWJfv7kudzx5N0b7WBmJ919btTnSMWgFfrLD8dZCIxVEli0c2XZ58u7rrtalUJCPQjsaIWqrX1jNRor+kIp+3x51/3F2rmLGpndcfNV+v7ibp1Y2EVQnxDk2NEKVWfcTTQay0qlFD1vtuzzDbtu6u0ZUD8CO5LWC5Z5K0FFZ9x1L4ZKylyc/ePfmdWXTi6PfKGU3WQ0zouKzUqTg8VTJCtrgbBfjIXAsouWsxsBdFRAzbtu3ovh4N4dkoq/qIou4uYFf14KaSi6eEpgR7LyKlmk5wNm08Fl2JiymKTHFncX+mxW8Ow1IRtUtppnWFVQ/336NzdJo18uMYP7JL5sigZ2UjFIVl5+2aSgJYplAkTZKpoyi7NZufH9R04HGUfe53spo17Qztqp+k8PPK5nBiaAsbtH0qp4OKpikKwmKlnKlhmWuXeIxdlQ/wZ5n99kNrLZ2WBQ74m5gzWVA9BTRWBHspo4XLpsgMgaU5ZQNeOh/g3yrpMXtPttsqxjjePuYKVV8XCkYpCsJo7MKxsg+sc0LP9fNVXUnx562fSULp26RKtPnhv73yDv33LYc0jDc+xZL5em8t5Fy0gnFYEdSau7JnucANEbU16lSdVvFIPXXV07p+mpTbrj5qsq/Vvk/VsOPkNvAbV/gXruVb9RurKnzrx3Sgegp4jAjolWJUDU9Y1iVHoo5P2KPsOwF2z/EYKD6lpk5QD04Sh3xMRLrWxu28I9uRuypqc2JVV2OGqvgVSu5BPDUe4IFJTaFvy89FBWBUvsssMiRwiS924eVTFAYspWsMSsBBl1b/LecRDYgQYVOcWp/6i7XofG/qPvBsWcEQ+7N22C46mUijGzQ5L+SNLTkv5X0p+5+2qIgQFdU6ZqpGgFS+wZcd7iMwE9rqo59nsl3ebu583sbyTdJukvqg8L6J5h1S5FgmDTlSBFFpVTHBMqBnZ3/1rfj1+X9M5qwwG6K8RuyaYWekN8u4g5pkkXMsf+fklfDXg9oFNineI0jhR7saQ4plSNDOxmdp+ZfSfjPzf1feZ2Secl3TnkOvvMbMnMllZWVsKMHmiRJnrfhJJiL5YUx5SqkakYd7922J+b2fsk3SDpGh+y28ndD0s6LK1vUCo5TqD12rRbMsVeLCmOKVVVq2Ku1/pi6Vvc/ckwQwK6K7XNUHlS7MWS4phSVbUq5lOSfk3Svbbe2vPr7v7nlUcFIKoUv12kOKZU0SsGAFqiaK8Ydp4CQMcQ2AGgYwjsANAxBHYA6BgCOwB0DIEdADqGE5QAtB5dHy9EYAfQanR9vBiBHZhAXZrhlulz36XnHobADkyYrs1wi3Z9DPncqb8gWDwFJkzX+poX7XMf6rl7L4jl1TW5nn9BZJ1fGwuBHZgwXetrXrTPfajnbsOLkcAOTJimTnI6empZ84vHtW3hHs0vHq9tRrtn56wO7t2h2ZlpmaSZ6SldOnWJ9h85fcF9Qz13G16MBHZgwjRxklPT6Yo9O2d1YmGX7rj5Kj11/ln9/MlzF9031HOP84Jo6iXXQ2AHJszgDHd2ZloH9+4IuvgXK10xqkImxHOXfUHEyMlTFQNMoLpPcgqZrihTgTLqviGeu+yBH2XKMUMhsAMILtT5pGVLFJs6F7XMCyJGTp5UDIDgQuWzy6Z0mlg/KKupxep+BHYAwYXKZ5ed7TaxflBWjJcNqRgAtehPV/Ty5PuPnC61U3Oc1Erd6wdlxTiEm8AOoFZVtvIfuG77BX9Xip9aGUfTLxtSMQBqVaX0McXUShswYwdQq6pVIamlVtogSGA3s49KOiRps7v/NMQ1AXRDUyWIZfXXx79sekpm0uqT55Ls1lhW5VSMmV0h6W2Sflh9OAC6JsUSxMHdoKtr5zLbELRViBz7HZI+JskDXAtAx6SYJ8/K+/dLrVtjWZVSMWZ2o6Rld/+WmQUaEoCuSS1PXiS/n1K3xrJGBnYzu0/SKzP+6HZJH5f09iI3MrN9kvZJ0pVXXlliiAC6JIXTh/Ly/oOfaauRgd3dr836vZntkLRNUm+2vkXSN83sanf/ccZ1Dks6LElzc3OkbYAJlMqxfFn18f3KrAGk8KIaNHaO3d0fdPdXuPtWd98q6aykN2YFdQCQ0jl9KOtwjpe/aKr0GkCqx+RRxw6gMSmdPhQi75/3orr1yGkdOnYm2uw92M7TjZk7NewAcsXodFinYS+kmLN3WgoAaMywmvamj48LYdQLKVbZJIEdQGPyatolJZmrHiXrRTUoRpqJHDuARmXltucXjzd+fFwI/S1588onY6SZmLEDiC6lRdWy9uyc1YmFXfrkzVcl0zqBwA4gui4sqqbUOsHcm98rNDc350tLS43fF0CaBjcuSZJpvQHVbMObflLccNRjZifdfW7U58ixA4huMFfdC+pSs7tTU9kZWxWpGABJ6OWqZ2emL2oV21TZYCo7Y6sisANISsyF1DYv4vYjsANISsyF1C4s4koEdgCJiXniUoqnPY2DxVMASelfSG26MiXmvUOi3BEAWqJouSOpGADoGFIxAKJJeTNQmxHYAUTRlc1AKSIVAyCKrmwGShGBHUAUXdkMlCICO4AourIZKEUEdgBRdGUzUIpYPAUQRds2A7WpgofADiCarGPyUtS2Ch5SMQAwQtsqeCoHdjP7sJmdMbOHzOxvQwwKAFLStgqeSqkYM3urpJskvcHdnzKzV4QZFgCk4/KZaS1nBPFUK3iqztg/IGnR3Z+SJHf/SfUhAUA4R08ta37xuLYt3KP5xeM6emq59DXaVsFTNbC/RtLvmdkDZvafZvamEIMCgBB6i57Lq2tyPb/oWTa479k5q4N7d2h2Zlqm9QO2D+7dkeTCqVQgFWNm90l6ZcYf3b7x918u6c2S3iTpC2b2as/oBWxm+yTtk6Qrr7yyypgBoJBhi55lg3LZCp6Y5ZEjA7u7X5v3Z2b2AUl3bwTy/zazZyVdJmkl4zqHJR2W1vuxjz1iACgo1qJn7PLIqqmYo5J2SZKZvUbSCyX9tOqgACCEqm0Lxs3Pxy6PrBrYPyvp1Wb2HUl3SXpfVhoGAGKosuhZJT8fuzyyUrmjuz8t6T2BxgIAQVVpW1AlPx+7PJKWAgA6bdy2BVVm3Qeu235Bjl1qtjySlgIAkKFKfj52eSQzdgDIUHXWHbPBGYEdADK0ra1wPwI7AORoS1vhQeTYAaBjCOwA0DEEdgDoGAI7AHQMgR0AOsZitHYxsxVJP2j8xtVdpslrcjZpzzxpzyvxzG3yKnffPOpDUQJ7W5nZkrvPxR5HkybtmSfteSWeuYtIxQBAxxDYAaBjCOzlHI49gAgm7Zkn7XklnrlzyLEDQMcwYweAjiGwj8nMPmpmbmaXxR5LnczskJk9YmbfNrN/NrOZ2GOqi5ldb2ZnzOxRM1uIPZ66mdkVZna/mT1sZg+Z2S2xx9QEM9tkZqfM7Muxx1IXAvsYzOwKSW+T9MPYY2nAvZJ+293fIOl/JN0WeTy1MLNNkj4t6R2SXi/pT8zs9XFHVbvzkj7i7q+T9GZJH5yAZ5akWyQ9HHsQdSKwj+cOSR+T1PkFCnf/mruf3/jx65K2xBxPja6W9Ki7f2/jLN+7JN0UeUy1cvcfufs3N/77r7Qe7NrXo7YEM9siabekz8QeS50I7CWZ2Y2Slt39W7HHEsH7JX019iBqMivp8b6fz6rjQa6fmW2VtFPSA3FHUrtPan1S9mzsgdSJgzYymNl9kl6Z8Ue3S/q4pLc3O6J6DXted/+Xjc/crvWv7nc2ObYGWcbvOv+NTJLM7MWSviTpVnf/Zezx1MXMbpD0E3c/aWZ/EHs8dSKwZ3D3a7N+b2Y7JG2T9C0zk9bTEt80s6vd/ccNDjGovOftMbP3SbpB0jXe3frYs5Ku6Pt5i6QnIo2lMWY2pfWgfqe73x17PDWbl3Sjmf2hpEslvdTMPu/u74k8ruCoY6/AzL4vac7d29hMqBAzu17SJyS9xd1XYo+nLmb2Aq0vDl8jaVnSNyS9290fijqwGtn67ORzkn7m7rfGHk+TNmbsH3X3G2KPpQ7k2DHKpyS9RNK9ZnbazP4h9oDqsLFA/CFJx7S+iPiFLgf1DfOS3itp18b/tqc3ZrNoOWbsANAxzNgBoGMI7ADQMQR2AOgYAjsAdAyBHQA6hsAOAB1DYAeAjiGwA0DH/D8RMkdliKjW8wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.linspace(-5, 5, 100), forward_sample['y'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_prob = model.make_log_prob_function()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: id=5237, shape=(), dtype=float32, numpy=-5309.3237>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_prob(forward_sample['sigma'], \n",
    "     forward_sample['intercept'],\n",
    "     forward_sample['x_coeff'],\n",
    "     forward_sample['y'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Centered-Eight "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'mu': <tf.Tensor: id=4177, shape=(), dtype=float32, numpy=4.7133865>,\n",
       " 'tau': <tf.Tensor: id=4205, shape=(), dtype=float32, numpy=0.12147237>,\n",
       " 'theta': <tf.Tensor: id=4227, shape=(8,), dtype=float32, numpy=\n",
       " array([0.70645565, 0.6838048 , 0.59054464, 0.611758  , 0.81999266,\n",
       "        0.7988661 , 0.71181107, 0.65780187], dtype=float32)>,\n",
       " 'y': <tf.Tensor: id=4249, shape=(8,), dtype=float32, numpy=\n",
       " array([ 2.2448833 ,  0.22619787,  0.31986347,  0.39096758,  1.4634166 ,\n",
       "        -0.25934392,  0.18226412,  0.62499356], dtype=float32)>}"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "schools = np.array([28.,  8., -3.,  7., -1.,  1., 18., 12.], dtype='float32')\n",
    "sigma = np.array([15., 10., 16., 11.,  9., 11., 10., 18.], dtype='float32')\n",
    "n_points = 8\n",
    "\n",
    "@pm.model(auto_name=True)\n",
    "def centered_eight(n_points = n_points):    \n",
    "    # Define priors\n",
    "    mu = pm.Normal(mu=0, sigma=5)\n",
    "    tau = pm.HalfCauchy(beta=.1, bijector=bijectors.Identity())\n",
    "    theta = pm.Normal(mu=tf.fill([8], mu), \n",
    "                      sigma=tf.fill([8], tau))\n",
    "    \n",
    "    # Define likelihood\n",
    "    y = pm.Normal(mu=theta, sigma=theta)\n",
    "    \n",
    "model = centered_eight.configure()\n",
    "\n",
    "forward_sample = model.forward_sample()\n",
    "\n",
    "forward_sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: id=4496, shape=(), dtype=float32, numpy=-5692.7417>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logp_func = model.make_log_prob_function()\n",
    "logp_func(forward_sample['mu'], forward_sample['tau'],\n",
    "          forward_sample['theta'], schools)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
